Aerodynamic Stability of Satellites in Elliptic Low Earth Orbits 



Matthew Bailey, Brian Ladson 

Stefan C. Mancas, Bogdan Udrea, Uchenna Umeadi 
Embry- Riddle Aeronautical University, 
Daytona Beach, FL. 32114-3900, U.S.A. 

Topical observations of the thermosphere at altitudes below 200 km are of great benefit in ad- 
vancing the understanding of the global distribution of mass, composition, and dynamical responses 
to geomagnetic forcing, and momentum transfer via waves. The perceived risks associated with 
such low altitude and short duration orbits has prohibited the launch of Discovery-class missions. 
Miniaturization of instruments such as mass spectrometers and advances in the nano-satellite tech- 
nology, associated with relatively low cost of nano-satellite manufacturing and operation, open an 
avenue for performing low altitude missions. The time dependent coefficients of a second order 
non-homogeneous ODE which describes the motion have a double periodic shape. Hence, they will 
be approximated using Jacobi elliptic functions. Through a change of variables the original ODE 
will be converted into Hill's ODE for stability analysis using Floquet theory. 

We are interested in how changes in the coefficients of the ODE affect the stability of the solution. 
The expected result will be an allowable range of parameters for which the motion is dynamically 
stable. A possible extension of the application is a computational tool for the rapid evaluation of 
the stability of entry or re-entry vehicles in rarefied flow regimes and of satellites flying in relatively 
low orbits. 

I. INTRODUCTION 

Non-dimensional aerodynamic coefficients are a measure of performance of a vehicle moving through a 
fluid at different speeds. The non-dimensional aerodynamic coefficients salient for this proposal are the 
pitch damping {Cma + C'm,) and pitch stiffness Cm„ coefficients. The pitch coefficients determine the pitch 
stability of the vehicle. 

In general, the pitch damping and stiffness coefficients are functions of the geometry of the vehicle, angle 
of attack, the rate of the angle of attack, and the pitch rate. For a vehicle operating in flow regimes with a 
Knudsen number larger than unity particle-surface interactions have a large influence on the pitch coefficients. 
These interactions are defined by energy accommodation coefficients. Examples of vehicles operating in this 
regime are relatively low flying spacecrafts and re-entry vehicles. 

In the preliminary stages of a low flying or re-entry mission study, a large level of uncertainty exists in 
the values used for the non-dimensional aerodynamic coefficients. The reasons for the high level uncertainty 
are multiple. The most important reason is the fact that little it is known about the energy accommodation 
coefficients. The second most important uncertainty is the composition of the atmosphere throughout the 
time interval of interest. 

The purpose of this paper is to analyze the pitch stability of a vehicle flying in the rarefied flow regime 
and determine ranges of pitch damping and pitch stiffness coefficients for which the pitch motion is stable. 
For this purpose we use analytical methods based on Floquet theory j^, and numerical methods developed 
in MATLAB fS^ to solve the following second order linear ODE 

a{t) - {Mg{t) + M^{t))a{t) - M^{t)a{t) = 9{t), (f) 

where a{t) is the angle of attack, M is the pitching moment, and 9 is the equivalent of the pitch angle 
and they are all functions of time. The following notation has been used: Mq = 77^, = j^^, 

Mq, = J , J being the moment of inertia. The relationship between the M coefficients of the ODE and 
the non-dimensional pitching moment coefficients is given in Eq. Q. The coordinate systems employed in 
the definition of the angles are presented in Fig. [f] 

Due to their observed double periodic shape, the periodic coefficients Mq{t) + Ma{t) (the damping term), 
and Ma{t) (the stiffness), will be approximated using Jacobi elliptic functions sn{(^ , m) , cn((^ , m) / dn{(^ , m) , 
where independent variable C is a function of time, and m is the modulus of the Jacobi functions used [2] . Of 
interest is how changes in the damping and stiffness coefficients will affect the stability of the pitch motion. 
The expected result will be an allowable range of parameters for which the motion is dynamically stable. 
One possible application of the technique is a tool for the rapid study of the stability of entry or re-entry 
vehicles. 




FIG. 1; Coordinate systems employed in the derivation of the equation of describing the pitching motion of the 
spacecraft. 



II. PROCEDURE 



A. Orbit 



This section refers specifically to the analysis of CubeSat mission. One of the mission's requirements is 
that the Dipping Thermospheric Explorer (DiPTE) CubeSat shall operate in an elliptical orbit of 700 km 
apogee and 200 km perigee. The other mission requirement relevant to this proposal is that the pointing 
knowledge of the payload sensitive axis in the ram direction shall be within 0.1° (Scr) full cone and the 
pointing control shall be 2.0° (3ct) full cone in the ram direction. 

The orbit requirement is derived from the need to perform measurements of density and temperature and of 
wind direction and magnitude of relevance for the propagation of various types of waves in the thermosphere. 
The attitude knowledge and control requirement is derived from the accuracy constraints of the payload. 



B. Approximation of Aerodynamic Moments 



It has been noticed numerically that the time dependent coefficients of the second order non-homogeneous 
ODE which describes the motion of the satellite have a double periodic shape. 

The forcing function 9{t) of Eq. ([!]) is obtained from the equations for the components of the velocity 
vector expressed in the satellite body frame and the satellite normal frame, and is due to the rotating frame 
in which the angle of attack was defined. This function is known and is also periodic and obtained from 
taking twice the time derivative of 

o/ X / cos + e cos 21/ \ 

eM = arctan -. — — , 2) 

where v is the true anomaly, and e is the eccentricity of the orbit given by e = ^^2£2_^j?££n ^ Here, the 
distance from focus to apoapsis Rapo = Radsarth + 700, and the the distance from focus to periapsis Rperi = 
Radsarth + 200, with the radius of the Earth being the equatorial Earth radius Rsarth — 6378.1363 -ftTm. 



In order to accomplish the overah goal of analyzing the pitch stability of this space vehicle flying in the 
rarefied flow regime and ultimately determine the ranges of pitch damping and pitch stiffness coefficients for 
which the pitch motion is stable, we used a preliminary method of an improved approximation of the density 
data using a least squares approximation of Jacobian elliptic functions [21 \5\ . 

Let Tji = y{ti) be the set of {n density p{t), or any other periodic function) data points that we are trying 
to model using a general least squares model. Then 



3=0 



Zj + e, 



(3) 



where Zj = sn{Q^ Kj) are m + 1 base elliptic functions, Q = j^ot are the harmonics, Kj the modulus of the 
elliptic sine function used, and e is the error that we want to minimize. 

Note that instead of sn, one could use cn,dn or any other combination (that works). It is important to 
have analytical expressions for the moments Mq{t), AIa{t), and Ma{t) since these coefficients determine the 
pitch stability of the vehicle, and they are needed to compute the angle of attack a{t) by solving Eq. ([T]). 

Using this method and having analytical expressions for p{t) in terms of elliptic functions, the damping 
and stiffness moments are obtained via 



1 



Mait) - -^j;^Qit)SreflrefC^m.i 
Ma{t) — —Q{t)SreflrefCm^, 



(4) 



where Q{t) 



is the dynamic pressure, p(t) is the mass density of the air, v{t) is the speed, Iref is 



the reference length (0.1 m), and Sref is the reference surface area (0.01 m ), and , Cm^, are the 
pitch damping and stiffness coefficients ( with Cm, + Cm^ = ~10, and Cm^ — —1.635 respectively). 
Using this notation, Eq. ([s]) can be written in matrix form as 



Y = ZA + E, 



(5) 



where the bases matrix Z is given by 



Z = 



sn{uJoti,Ki) 
sn{uJot2, Ki) 



1 sn{uJotn,Ki) ■■■ sn{muJotn,Km). 



(6) 



Therefore, the sum of squares of residuals is formulated as follows 



n 

E(. 



(7) 



By taking m + 1 partial derivatives with respect to Uj, the normal equation is A = {Z*^ Z)~^ Z'^Y , where 
A = (flQ ai ... a,„)* are the unknown coefficients of the base functions that we are trying to find, and 
y — (j/i 2/2 ■•■ UnY are raw the density data points. 

In order to accomplish a least squares model of the periodic density function, first we determine the period 
of the data which is T = 5615.2 sec. This value was determined from the original Pitch Dynamics MATLAB 
code. Using this value for period, ujq is found from ujq = 2tt/T. 

Let the Jacobian elliptic integral [2] be defined as 



then X ~ sn(u, k). When x — 1^ the complete elHptic integral becomes 
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K^r (9) 

which has the quarter period K — T/A. Because the above integral that defines the full elliptic sine function 
cannot be solved analytically, MATLAB was used to systematically find the modulus that matched with the 
raw data's period of T = 5615.2 sec. The modulus found using Eq. ^ was k — 0.4405. In this procedure 
we assumed that all base functions have the same modulus k = k"^ . 

The data used for approximating the density as an elliptic sine approximation comes from a mass density 
profile text file developed by Rick Doe, SRI International [Ij. The MATLAB code reads in the altitude and 
density columns for the file in order to conduct proper calculations using the appropriate data. However, 
before calculations begin, the altitudes of interest are extracted with their corresponding density values. Of 
the original 486 data points from the text file, only 251 points arc useful. These 251 points acquire half the 
period, which means there are a total of 502 data points to fill one orbit. The true anomaly values range 
from to 359 degrees, but we wanted to evaluate our motion starting from the apogee. Therefore, tt must 
be added to each true anomaly value (after being converted into radians) in our range. Based on our data 
points, a step size is created accordingly to fit the data. This step size is determined by taking the largest 
value of the range (359 degrees) and dividing it by one less the total number of data points (501). The 
altitude at the perigee and apogee is 200 km and 700 km, respectively. The radius of the Earth at these 
altitudes can be determined, which can lead to values for the semi-major axis, eccentricity, and orbit period. 
The orbital position and speed provide calculations using Kepler's formulas, and this will be explained later. 

In order to begin the elliptic sine approximation, the data points were sorted accordingly based on the 
fact that motion is starting from the apogee. The Z matrix is developed via Eq. (|6]). As discussed above, 
the A matrix is [Z^Z)^^Z^Y . Eq. ([5]) develops the matrix for the density approximation using elliptic 
sine functions. The sum of the squares residual that was calculated using the elliptic functions method was 
Sr = 7.921010"^^ which is a far better method of approximation of the data than the use of the cubic spline 
approximation method that had an error Sr = 2.683910~^. 

By consulting the comparison in Fig. [2] it is notable that the percent difference between the errors 
calculated for the two methods was 199.988%. Regardless that the percent difference error is so large, 
the errors of these approximations are relatively small in comparison to the raw data on the logarithmic 
scale. 

Since we have analytical expression for the density, by using Eq. Q, we obtained the moments Mq{t), 
Ma{t), and Ma{t). The comparison between the cubic spline approximation and the elliptic sine approxi- 
mation of the aerodynamic moments are shown in Fig. [3j 

C. Pitch Dynamics 

The 700 km x 200 km elliptic orbit makes the pitch dynamics an interesting dynamic problem due to the 
fact that the density varies by three order of magnitudes between the apogee and the perigee. The density 
at the apogee is so low that the dynamics of the pitch motion is very close to that of a double integrator. 
At the perigee the density is sufficient to provide enough pitch stiffness and some pitch damping. 

Numerical integration of the homogeneous ODE ^ using MATLAB, describing the pitch motion shows 
the expected behavior. The results of the integration are presented in Fig.|4] The moments used initially are 
based on the cubic spline approximation. The initial conditions used are pitch angle of ao = 2° and pitch 
rate of do = 0. The pitch settles in an oscillatory motion with oscillations at the perigee with an amplitude 
of 0.35° well within the 2° full cone requirement. It is interesting to take a look at the pitch motion in the 
phase plane, i.e., pitch angle vs. pitch rate, shown in Fig. [5] The phase plane diagram is similar to that of a 
limit cycle motion. Since the motion of the pitch angle and pitch rate indicates a quasiperiodic motion, in 
the phase plane the motion is limited to an attractor. 

D. Pitch Stability 

The purpose of this paper is to analyze the pitch stability of the CubeSat flying in the rarefied flow regime 
and to determine its ranges of pitch damping and pitch stiffness coefficients for which the pitch motion is 
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FIG. 2: Density raw data (top), approximation using splines (middle), approximation using elliptic functions (bottom). 



stable based on elliptic functions. In order to do so, we performed an analysis using analytical methods 
delineated in Floquet theory [3|, together with numerical methods developed in MATLAB. We will now 
discuss the methods and algorithms utilized in the construction of operational MATLAB code in order to 
solve the nonhomogeneous (ODE) as denoted in Eq. 0. 

First, after many different tries, we have established that the most appropriate ODE solver script for our 
problem is provided was the odell3 solver. This is primarily used for solving non-stiff differential equations 
by the means of variable order method. This solver integrates the system of differential equations over the 
time period given by tspam and evaluates the solution starting with the initial conditions given by yo- The 
output is the column vector yoDE with corresponding time vector tpDE- 

To find the right-hand side of Eq. ([T]) , we differentiated twice Eq. ^ and we obtained 



0{i>, v) = arctan I — 



1 + 2e^ + 3e cos v 



1 



2e cos V 



(10) 



B{v, z>, v) 



1 + + 3e cos v . 



1 



e(e2 



1) sin V 



2e cos V 



(1 



2e cos v)- 



(11) 



It is important to note that the true anomaly v{t), along with its first two time derivatives, are also 
functions of time. Their expressions will also be approximated using elliptic sine functions as follows. Since 
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FIG. 3: Aerodynamic moments, approximation using splines (left), approximation using elliptic functions (right). 



via Eq. ^ 



then by differentiating, 



m 



(12) 



(13) 



and 



i) = ujq ''^^ aj sn{jijjQt)[2k'^ sn^ (joJot) — (1 + fc^)]. 



(14) 



The values of cos v, and sin v necessary in obtaining were calculated using the Keplerian equations that 
govern elliptically orbiting bodies. From the equation of radius R 



R = 



a(l -e^) 
1 + e cos v ' 



(15) 



where the semimajor axis, we find cos v and subsequently sin ly, and we substitute them 



in Eq. (11). 251 data points data points delineating the evolution the satellite altitude from its apogee 



to its perigee were concatenated in order to represent the altitude evolution for a full orbit and then were 
populated into a vector that served as a raw data for our elliptic sine approximation. Furthermore, the speed 
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FIG. 4: Time history of the pitch angle and pitch rate for 50 orbits (top) and a zoom in for the last five orbits 
(bottom). 



v{t) was determined via Eq. (16) 



where /i = 3.986004415 10^ ^™ represents the Earth's gravitation parameter. 

Using the analytical expression for i>{t) in terms of the elliptic sine function, Eq. and the MATLAB 
odellS solver, we solved numerically both the homogenous and nonhomogenous ODE Eq. 0. The resuhs 
are graphically plotted and analyzed below. Fig. [6] shows the homogeneous solutions of a{t) and a{t) over a 
full 10 periods (orbits). The figure actually shows an overlay of the solution as determined using cubic spline 
approximations and elliptic sine approximations. It is important to note that because these solutions are 
homogeneous, (the right hand side of the ordinary differential equation is null) they do not take into account 
the true anomaly approximations which solely appear in the forcing function 9. It is evident, as previously 
noted, that the newly found results from the elliptic sine approximations closely simulate the cubic spline 



Pitch Angle and Pitch Rate 8 




-0.02 -0.015 -0.01 -0.005 0.005 0.01 0.015 0.02 



d c: / dt (deg/s) 

FIG. 5: Phase plot diagram of the pitch dynamics. 



approximation. This is a testament to the accuracy of the results and the rehabihty of the approximation 
method as employed. Fig. [6] also testifies to the accuracy of the results and the consistency of the findings. 

Fig. [t] shows the phase plot overlay of the same information (the homogeneous solution from the cubic 
spline and elliptic sine approximations). The behavior of both approximated solutions were similar. The 
shape and predicted envelope of the pitch angle and pitch rate are comparable. The solutions produced by 
the two methods of interest begin to diverge; however when the non-homogeneous solutions are explored, the 
subsequent figures show that these solutions and the variations are more evident. At first glance, it is readily 
noticeable that both solutions exhibit the same general characteristics, namely a double periodic behavior 
with beats, as seen in Fig. [8j The simulation of longitudinal dynamics shows a damped oscillatory behavior 
with these beats that line up without a phase shift that would otherwise suggest inaccuracy in the results. 
Although the amplitude of the elliptic sine results is consistent with the previous findings with respect to 
the comparable cubic spline results in that they are larger. These larger amplitudes can be attributed to the 
more direct approximation that is done through the elliptic sine function that yields results that arc more 
accurate yet less controlled. 

Looking now at the phase plot overlay of the same results, we can develop a similar conclusion. Though 
the general attitude (behavioral shape) of the competing solutions is similar, there is a distinct issue of 
difference when it comes to the magnitude of the results as displayed in Fig. [9j The elliptic sine results yield 
values that are on a scale of almost a factor of 2 both for the pitch angle and the pitch rate. This indicate 
an error in the calculations, such as derivations for j> and i), or may just be the reality of what the elliptic 
sine approximations yield. The final conclusion is yet to be made as the researchers continue to investigate 
the results and possible alterations to the method used to solve the ODE. 

The results indicate a strong correlation between the cubic spline and elliptic sine approximations, yet do 
suggest that there may be discrepancies in the method used for finding and approximating the true anomaly 
values that are employed in the resolution of the ODE via elliptic sine methods. 

Further investigation will be considered on the use of other Jacobian periodic elliptic functions in order 
explore the aerodynamic longitudinal pitch stability of the satellite. Floquet theory is the next step of the 
comprehensive process in developing a full analytical and numerical evaluation of the stability regimes of the 
CubeSat. 



E. Dynamical Systems Analysis Using Floquet Theory 

The parameters associated with the dynamics of the motion of the of a CubeSat class mission flying in a 
700 km X 200 km orbit have been computed with a direct simulation Monte Carlo code or extrapolated from 
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FIG. 6: Pitch angle and pitch rate overlay for the homogeneous ODE. 
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FIG. 7: Phase plot overlay of the 2 methods solving the homogeneous ODE. 
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FIG. 8: Pitch angle and pitch rate overlay for the nonhomogeneous ODE. 
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FIG. 9: Phase plot overlay of the 2 methods solving the nonhomogeneous ODE. 



existing data. In this project we also determined the stabiHty regions of Eq. ([T]) describing the one 
of freedom attitude dynamics in low altitude elliptic orbits using Floquet therory. 
Once the damping and stiffness coefficients are known, by the change of variables 



(17) 



the homogeneous Eq. ([T]) will be converted into Hill's equation, Eq.(18). This equation was introduced by 
George W. Hill in his studies of the motion of the Moon. Roughly speaking the motion of the Moon can be 
viewed as a harmonic oscillator in a periodic gravitational field. Since the analytic solutions of Eq. ([T]) are 
not known, to analyse the stability, we will be using Floquet theory. Here, b{t) is also a periodic function 
such that b{t) = b{t + T) for some T which needs to be found; and it will act as an energy source for the 
system. 



m + h{t)p{t) - 



(18) 



In order to solve Eq. (18), h{t) must be determined first. Using Eq. Q, set oi = —[Mq{t) + Ma{t)] 
and a2 — —Ma{t) in order to apply the transformation given by (17). The ODE now possesses the form 



a{t) + aia{t) + a2a{t) = 9{t). After developing the appropriate derivatives of Eq. (17), substituting into the 
transformed ODE, and completing some algebraic manipulation, the transformed ODE yields the expression 



provided by (18) such that 



h{t) = -A'Ut) + ~[M,it) + AUit)] - ^[M,{t) + A'U{t)]' 



(19) 



9{t) will transform to 7(i), however Hill's ODE solves for a homogenous differential equation. Therefore 
7(i) does not need to be determined in this case. To continue with proving the stability of the satellite, b(t) 
must be calculated. According to (|4]), 



q 7 ^ 

Oref'ref 

4J 



+ a„ jt'(i)p(<) 



(20) 



where ■^[v{t)p{t)] — v{t)p{t) + v{t)p{t). Using (20) and the expression for AIa{t) from and letting an 
arbitrary constant A — '^''^^'''^^ 



4J 



, then 
A 



b{t) = ^2A^p{t)[v{t)f + ^(C™, + c„.jmp{t) + v{t)m] 

f'ref ^ 



(C„.,+C„J^[i;(t)]2[p(i)]^ (21) 



Furthermore, p(t) is obtained via ([s]) and by using the identities, ^sn(u) — cn{u)dn{u), and cn{u) = 
\/l — s'n?{u), dn{u) = yjl — k'^sn?{u), i.e., see also Eq. (13), with v replaced by p. 
Using linear stability of dynamical systems, we will transform Eq. ( 18 ) into a first order equivalent system 










1 ' 
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^b{t) 







. ^2 . 



(22) 



Applying linear systems theory, the stability of the zero solution of the linear periodic system can be analyzed. 
From the scalar second order linear differential equations, if il + p{t)u + q{t)u = 0, and the Wronskian of the 
two solutions ui and U2 being defined as 



W{t) det 



Ui{t) U2{t) 
Ui{t) U2{t) 



(23) 



then by Liouville's Lemma, 



W{t) = iy(0)e--/'oP(^> 



(24) 



Let the characteristic multiphers for Hill's Equation be denoted by pi and p2 and note they are roots of til 
characteristic equation 



p2 - (tr$(r))p + det$(r) = 0. (25) 
For notaional convenience let us set 2(/) — tr<i>(T), to obtain the equivalent characteristic equation 

p2 - 2(/)p + 1 = (26) 

whose solutions are given by 

pi,2 = 0± (27) 
Even though the solutions u are not known explicitly, we have that the Floquet multipliers pi^2 satisfy Eq. 



(261 where cj) — ^{ui{T) +U2{T)). The characteristic Floquet exponents /^i.2 are given by pi^2 — e^^-^"^, and 
consequently using Eq. ( 26 ) 



Ail + = 

coshpiT^cf) (28) 
Although is not known explicitly, it is useful to characterize the properties of pi,2, or /ii.2 in terms of cf). 



Hence, by Eq. (27) we have the following cases: 

(i) if (/) > 1, then pi 2 are distinct positive real numbers such that pip2 = 1- Thus, we may assume that 
< pi < 1 < P2, with pi — l/p2, and there us a real number ^ > (a characteristic exponent) that 
e^"^ = p2 and e"^"^ = pi. Then, there is a fundamental solutions set of the form e~'^*pi(i), and e'^*p2(t) 
where the real functions pi.2(t) are T periodic. In this case the zero solution is unstable; 

(ii) if (/> < —1, then pi 2 are both real and both negative. Also, since piP2 — 1 then we may assume that 
pi < —1 < p2 < 0, with pi — 1//32- Thus, there us a real number /i > (a characteristic exponent) 
such that e^^"^ = pf and e^^^"^ — P2- in the case (i), there is a fundamental solutions set of the 
form e''*gi(i), and e^^*g2(0 where the real functions qi,2{t) are 2T periodic. Again, the zero solution 
is unstable; 

(iii) if —1 < < 1, then pi_2 are complex conjugates with nonzero imaginary parts. Since pipb = 1, we 
have \pi \ — 1, and therefore both characteristic multipliers lie on the unit circle in the complex plane. 
Because pi^2 have nonzero imaginary parts, one of this characteristic multipliers, say pi, lies in the 

upper half plane. Thus, there is a real number with < OT < ir and e*^"^ = Pi- In fact, there is a 
solution of the form e*^*(r(t) + zs(t)), where r{t), s{t) are both T periodic functions. Hence, there is a 
fundamental solutions set of the form r(t) cos{9t) — s{t) sm{6t), r{t) sin(0t) + s{t) cos{9t). In particular, 
the zero solution is stable but not asymptotically stable. Also, the solutions are periodic if and only if 
there are relatively prime integers m and n such that 2TTm — nOT . If such integers exist all solutions 
have period nT. If not, then the solutions are quasi-periodic. 

Certain curves of the form (p — ±\ separate parameters regimes where unbounded solutions exist, i.e., 
|(/)| > 1, from regions where all solutions are bounded , i.e., |(/)| < 1. We have just proved the following facts 
for Hill's equation, and the results are summarized in the following Lemma , [3], which we will act as the 
stability/instability criterion in determining the parameters' regimes for which the solution is stable/unstable. 



Lyapunov Lemma: // b{t)dt < 4, then all solutions of the Hill's equation (18) are bounded. In 
particular the trivial solution is stable. 

The development of the stability region for the satellite in low elliptic Earth orbit followed the Lyapunov 



Lemma is stated. Based on this, Eq. (21 ) was integrated over to T. Hence 



b{t) = -2A^ rp(t)[vit)rdt+^{Cm,+C„,J f ^[v{t)p{i)]dt (29) 

're/ Jo ^ JO 
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FIG. 10: Pitch- damping coefficients based on x/D for Mach 2 and 5. 



where p{t)[v{t)]^dt = Ii and [v[t)Y[p{t)]^dt = I2 are constants. Eq. (29) can be reduced to 



(30) 



For simplicity, let fci = -2y^/i, ^2 = |[i;(r)p(f )-t;(0)p(0)], and fcg = such that fciC^^ +fc2(C„,, + 

Cma) + k3{Cm + CmaY ^hc simplified expression for h{t). 

In order to successfully find the constants Ii and /2, some method of numerical integration is needed. 
Typically, integrations can be solved analytically if an expression exists and the integrand is integrable. In 
this case, a numerical method was used for the elliptic sine approximations for the density and speed which 
posed a problem to evaluating a conventional integral. Necessary numerical approximation methods for these 
integrals were considered, such as Boole's rule and the trapezoidal rule. 

Boole's rule is a variation on the Newton-Cotes' formula developed by George Boole [5|. This method of 
numerical integration approximates the integral of the type 



f{x)dx 



(31) 



This method evaluates the integral by using values of / over five equally spaced steps of xi^ X2 = xi + h, 



FIG. 11; StabiUty region based on the first constraint, Cm^ + Cm^ ~ —10. 




FIG. 12: Stability region surface. 



FIG. 15: Contour plot of stability region based on lower bound. 




FIG. 16: Contour plot of stability region based on elevation {Cma) with lower bound. 



X3 ~ xi + 2h, X4 ~ xi + 3ft., and x^ — xi + Ah. Eq. (31 1 now exists as 
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^5 2/i r 1 

f{x)dx = — 7 fix,) + 32/(x2) + 12f{x3) + 32/(.T4) + 7/(x5) + E (32) 
45 L J 

where the error term is £' = ~ 1*3.^5*7*9 h"^ f^^^ (c) for some number c € [xi, x^]. 

The trapezoidal rule is a method of numerical integration where a finite number of trapezoids are con- 
structed beneath the curve of the function of interest. The method involves finding the area of each of these 
trapezoids and then summing up these areas to find an approximation of the total area under the curve. 



This is accomplished by using Eq. ( 33 ) , which is a variation on the simple equation for the area of a 
trapezoid. The error is minimized in this case, by maximizing the number of points that are used to do the 
approximation. In this case, the researchers used 502 points, matching the number of data points from the 
original density and altitude data. 

f[x)dxM [b-a) (33) 

Like all numerical integration techniques, an error exists. An estimation for the local truncation error of a 
single application of the trapezoidal rule is — ]i^/"(C)(^ ~ '^)^ where ^ lies somewhere in the interval from a 
and b. This error term indicates higher order functions (with curvature) might consist of some error, while 
a linear function yields no error. The calculation of Ii and I2 utilized the trapezoidal rule for its simplicity. 

The first constraint, Cm, +Cm^ = —10, is a first order guess such that the pitch damping is similar. Pitch 
damping of flared projectiles has been investigated by Weinacht and co-authors [1]. The study assumed 

• long and slender axi-symmetric shape with a sharp nose 

• flight at sea level at speeds of 680 and 1700 m/s (Mach 2 and 5) 

• zero-spin coning motion, i.e., coning motion with null spin rate 

The Dipping Thermospheric Explorer (DipTE), however, is short and squat with a blunt nose and flics in 
a elliptic orbit at speeds of 7366 and 7950 m/s. With an aspect ratio (L/D) of 3, and ratio of the center 



of gravity to length (x/L) of 0.5, the ratio x/D is 1.5. From Fig. 10 it is determined that Cm^ ~ and 



Cm ~ —10. This is under the assumption that the center of mass is at the geometric center, i.e. the worst 



case scenario. This constraint is applied to Eq. ( 30 1 to determine an upper surface boundary under the 



Lyapunov Lemma. If the constraint is applied to the simplifled form of b(t), then kiCma ~ 10fc2 + lOO/ca < 



if fci > 0. Therefore, it is determined that C,„^ > —0.0283792. The representation of Fig. 11 displays the 
stability region of applying the first constraint. It can be seen that the blue plane represents a graphical 
model of this constraint, as the magenta plane represents the boundary plane as discussed. Therefore, any 
values above the magenta plane would yield stability as the values below are unstable. 

The second constraint assumes that and C™^ are free variables because the stability of the satellite 
depends on these coefficients. In this case, if their summation is equivalent to an arbitrary variable w, then 
/(ui) = fciCniQ +fc2W+fc3U'^ < 4 assuming /ca > 0. The necessary condition to have is A = ^2^— 4fc3(fciC„i„ — 
i) > which may be different from the first condition. This condition exists due to the quadratic nature of 
the inequality and then proves and provides the lower boundary such that C„iq < —0.0267126. Therefore, 
the stability region will need to fall below this value as well. 



Fig. 13 demonstrates a graphical representation of the comparison of b(t) to the criterion < 4 using MAT- 



LAB to ensure the bounded nature of Hill's equation. This figure shows a two-dimensional representation of 
the stability region of the satellite in terms on the pitch rate and pivot coefficients. This graphic is essen- 



tially a contour plot of the figure based on elevation from the three-dimensional model as in Fig. 12 The 
multiple colors are generated by MATLAB to show this difference in graphical altitude, where the lighter 
colors (dark blue, blue, light blue, yellow-green, orange etc.) represent the lower regions and the darker red 
colors represent the higher regions. Specifically in this case, we are permitted immediately to examine the 
stability region of the satellite and associate it with the dark red region. This region falls between the upper 
and lower bounded surfaces of the criterion associated with the transformation of Hill's equation. 

Fig. [14] shows a three-dimensional representation of the region in which the satellite is expected to be 
stable and by default the remaining region where the satellite is expected to be unstable. The figure is made 
up of three distinct surfaces. The curved surface, although appears flat on this scale, corresponds to the full 
domain and range of the stability coefficients Crn„ , , and Cm^ ■ This is the full range of possibilities for 



the combination of these coefficients as the satellite traces its orbit and as limited by the physical limitatio}i§ 
of the satellite as mentioned in the previous section. This figure shows an upper and lower curved boundary 
which displays the stability region lying between the surfaces. These boundaries were developed from the 
second constraint based on the quadratic expression f{w) used with the criteria based on the Lyapunov 
Lemma. As previously discussed, a Cm^ and Cm^ (-0.0267f 26 and -0.0283792 respectively) were 
determined. This range was used in order to develop the roots of the quadratic expression such that 



where A is a function of Cm^^ correspondingly making the roots wi^2 functions of Cm,- Since the roots 
do exists, i.e. f(w) < if and only if wi < w < W2- Likewise, Cm^ < W2 — y and Cm, > wi — y. 
Therefore the surface as seen in Fig.[T2]has a small portion representing the range and domain of the stability 
coefficients that would allow the satellite to remain stable between the boundary surfaces Cm, and Cm, ■ 
Since it is evident that the curved surface provided in Fig . |12| intersects the lower surface, a contour plot 
can be developed to provide the region of stability. In Fig. |15[ this plot gives elevation data for the lower 
boundary surface. The portion that the stability curve can be accounted for is the first dark blue region. The 
remaining colors are unstable regions due to falling below the lower bound. When combining the contour 



plots developed in Fig. 13 and Fig. 15 a final representation of the stability region based on elevation with 
the lower bound included can be shown by Fig. |16| The gray region shows where all the conditions are met 
for stability of the satellite. 



III. CONCLUSION 



In this paper we used analytical and numerical methods and determined the stability regions of the 
equation describing the one degree of freedom attitude dynamics in low altitude elliptic orbits. The time 
dependent coefficients of the second order non-homogeneous ODE which describes the motion had a double 
periodic shape. Hence, to approximate them we used a novel and powerful technique based on Jacobi elliptic 
functions using Jacobi elliptic sine function. Through a change of variable the original ODE which described 
the motion of the satellite was be converted into Hill's ODE suitable for stability analysis using Floquet 
theory. This allowed us to establish how changes in the coefficients of the ODE affect the stability of the 
solution via all the transformations. The expected result was be an allowable range of parameters for which 
the motion is dynamically stable or unstable. A possible extension of the application is a computational tool 
for the rapid evaluation of the stability of entry or re-entry vehicles in the rarefied flow regimes. 
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